SS3 Course ICES 2023
Mauricio Mardones
19 September, 2023
Background
This code contains model files and associated material for the 2023 ICES course on Stock Synthesis held in Copenhagen, Denmark.
Model setup
# install.packages("devtools")
# devtools::install_github("r4ss/r4ss", ref="development")
#
# install.packages("caTools")
# library("caTools")
#
# install.packages("r4ss")
library(r4ss)
library(here)
#remotes::install_github("PIFSCstockassessments/ss3diags")
library(ss3diags)
library(knitr)
library(tidyverse)dir1<-here("herring-aspm")# version 3.30.21
dir2<-here("Workshops" , "SS3_ICES_Course", "herring-just-2019-as-equilibrium")
dir3<-here("herring-scaa") #
dir4<-here("Workshops" , "SS3_ICES_Course", "herring-scaa-2dar") #
dir5<-here("herring-scaa-2dar-noTV")
dir6<-here("Workshops" , "SS3_ICES_Course", "predators") # v. 3.30.21Run models
#shell(cmd="ss_win") # run SS para windows
# aqui debo cambiar el setwd() para correr el modelo
#dir5 <- setwd("~/DOCAS/SA_Krill/s5")
# r4ss::run(dirvec = dir3, model = './ss_osx',
# skipfinished = FALSE)
# #system('./ss_osx')
#another way to run with new r4ss
r4ss::run(
dir = dir1,
exe = "ss_osx",
skipfinished = FALSE, # TRUE will skip running if Report.sso present
show_in_console = TRUE # change to true to watch the output go past
)
r4ss::run(
dir = dir1,
exe = "ss_osx",
skipfinished = FALSE, # TRUE will skip running if Report.sso present
show_in_console = TRUE # change to true to watch the output go past
)
r4ss::run(
dir = dir2,
exe = "ss_osx",
skipfinished = FALSE, # TRUE will skip running if Report.sso present
show_in_console = TRUE # change to true to watch the output go past
)
r4ss::run(
dir = dir3,
exe = "ss_osx",
skipfinished = FALSE, # TRUE will skip running if Report.sso present
show_in_console = TRUE # change to true to watch the output go past
)
r4ss::run(
dir = dir4,
exe = "ss_osx",
skipfinished = FALSE, # TRUE will skip running if Report.sso present
show_in_console = TRUE # change to true to watch the output go past
)
r4ss::run(
dir = dir5,
exe = "ss_osx",
skipfinished = FALSE, # TRUE will skip running if Report.sso present
show_in_console = TRUE # change to true to watch the output go past
)
r4ss::run(
dir = dir6,
exe = "ss_osx",
skipfinished = FALSE, # TRUE will skip running if Report.sso present
show_in_console = TRUE # change to true to watch the output go past
)Diagnosis
herring-aspm Model
Saco los outputs en html
Data disponible para este escenario. Espinel es la serie mas consistente del conjunto de datos.
SSplotData(base.model1, subplot = 2)## Warning: The `subplot` argument of `SSplotData()` is deprecated as of r4ss 1.45.1.
## ℹ Please use the `subplots` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Respecto a los valores y parametros biologicos modelados, los siguientes graficos identifican los estimadores puntuales del recurso
SSplotBiology(base.model,
subplots =1,
labels = c("Length (cm)", "Age (yr)",
"Maturity",
"Mean weight (kg) in last year",
"Spawning output",
"Length (cm, beginning of the year)",
"Natural mortality",
"Female weight (kg)",
"Female length (cm)",
"Fecundity",
"Default fecundity label",
"Year",
"Hermaphroditism transition rate",
"Fraction females by age at equilibrium"),
)aporte de las cohortes por año para las capturas.
SSplotCohortCatch(base.model1)
AJuste de tallas por flota
SSplotComps(base.model1, subplots = 1)Otros plots
SSplotDynamicB0(base.model1, uncertainty = F)## Plotting Dynamic B0
#SSplotSPR(base.model3)
SSplotPars(base.model1)## Plotting distributions for 4 estimated parameters (deviations not included).## $mcmc not found in input 'replist', changing input to 'showpost=FALSE'
Salida de las biomasas con las dos flotas
SSplotTimeseries(base.model1, subplot = 1)
SSplotTimeseries(base.model1,
subplot=1,
add = FALSE,
forecastplot = FALSE,
uncertainty = FALSE,
bioscale = 1,
minyr = -Inf,
maxyr = Inf,
plot = TRUE,
print = FALSE,
plotdir = "default",
verbose = TRUE,
btarg = "default",
minbthresh = "default",
xlab = "Year",
labels = NULL,
pwidth = 6.5,
pheight = 5,
punits = "in",
res = 300,
ptsize = 10,
cex.main = 1,
mainTitle = FALSE,
mar = NULL
)
herring-scaa Model
Saco los outputs en html
Data disponible para este escenario. Espinel es la serie mas consistente del conjunto de datos.
SSplotData(base.model3, subplot = 2)
SSplotRecdevs(base.model3)
Respecto a los valores y parametros biologicos modelados, los siguientes graficos identifican los estimadores puntuales del recurso
SSplotBiology(base.model3, subplots =1,
labels = c("Length (cm)",
"Age (yr)",
"Maturity",
"Mean weight (kg) in last year",
"Spawning output",
"Length (cm, beginning of the year)",
"Natural mortality",
"Female weight (kg)",
"Female length (cm)",
"Fecundity",
"Default fecundity label",
"Year",
"Hermaphroditism transition rate",
"Fraction females by age at equilibrium"),
)aporte de las cohortes por año para las capturas.
SSplotCohortCatch(base.model3)
AJuste de tallas por flota
SSplotComps(base.model3, subplots = 1)Otros plots
SSplotDynamicB0(base.model3, uncertainty = F)## Plotting Dynamic B0
#SSplotSPR(base.model3)
SSplotPars(base.model3)## Excluding 57 deviation parameters because input 'showdev' = FALSE## Plotting distributions for 14 estimated parameters (deviations not included).## $mcmc not found in input 'replist', changing input to 'showpost=FALSE'
Salida de las biomasas con las dos flotas
SSplotTimeseries(base.model3, subplot = 10)SSplotTimeseries(base.model3,
subplot=1,
add = FALSE,
forecastplot = FALSE,
uncertainty = FALSE,
bioscale = 1,
minyr = -Inf,
maxyr = Inf,
plot = TRUE,
print = FALSE,
plotdir = "default",
verbose = TRUE,
btarg = "default",
minbthresh = "default",
xlab = "Year",
labels = NULL,
pwidth = 6.5,
pheight = 5,
punits = "in",
res = 300,
ptsize = 10,
cex.main = 1,
mainTitle = FALSE,
mar = NULL
)herring-scaa-2dar-noTV Model
Saco los outputs en html
Data disponible para este escenario. Espinel es la serie mas consistente del conjunto de datos.
SSplotData(base.model5, subplot = 2)
Respecto a los valores y parametros biologicos modelados, los siguientes graficos identifican los estimadores puntuales del recurso
SSplotBiology(base.model5, subplots =1,
labels = c("Length (cm)",
"Age (yr)",
"Maturity",
"Mean weight (kg) in last year",
"Spawning output",
"Length (cm, beginning of the year)",
"Natural mortality",
"Female weight (kg)",
"Female length (cm)",
"Fecundity",
"Default fecundity label",
"Year",
"Hermaphroditism transition rate",
"Fraction females by age at equilibrium"),
)## Note: this model uses the empirical weight-at-age input.
## Plots of many quantities related to growth are skipped.## Warning in max(biology[["Len_mean"]]): no non-missing arguments to max;
## returning -Infaporte de las cohortes por año para las capturas.
SSplotCohortCatch(base.model5)
AJuste de tallas por flota
SSplotComps(base.model5, subplots = 1)Otros plots
SSplotDynamicB0(base.model5, uncertainty = F)## Plotting Dynamic B0
#SSplotSPR(base.model3)
SSplotPars(base.model5)## Excluding 69 deviation parameters because input 'showdev' = FALSE## Plotting distributions for 3 estimated parameters (deviations not included).## $mcmc not found in input 'replist', changing input to 'showpost=FALSE'
Salida de las biomasas con las dos flotas
SSplotTimeseries(base.model5, subplot = 10)SSplotTimeseries(base.model5,
subplot=1,
add = FALSE,
forecastplot = FALSE,
uncertainty = FALSE,
bioscale = 1,
minyr = -Inf,
maxyr = Inf,
plot = TRUE,
print = FALSE,
plotdir = "default",
verbose = TRUE,
btarg = "default",
minbthresh = "default",
xlab = "Year",
labels = NULL,
pwidth = 6.5,
pheight = 5,
punits = "in",
res = 300,
ptsize = 10,
cex.main = 1,
mainTitle = FALSE,
mar = NULL
)
Desembarques
Los desembarques utilizados para cada una de las flotas que generan remoción en el recurso reineta, a saber; espinel, enmalle e industrial (Figura 6).
Desembarques de merluza común por flota
Predator fleet with RW
Random Walk (RW) refers to a mathematical model that describes a stochastic process in which a variable changes randomly over time, without a clear trend or pattern.
Specifically, a random walk can be used as a Bayesian estimation technique to infer the posterior distribution of an unknown parameter. In this approach, it is assumed that the prior distribution of the parameter is a normal distribution with mean zero and a known variance, and that the parameter value at each time step follows a random walk process. Based on the observed data and the prior distribution, the posterior distribution of the parameter can be calculated using Bayesian inference.
Random walk is a useful tool for parameter estimation in dynamic models, as it allows for modeling uncertainty and variability in the parameter’s evolution over time. However, it is important to note that the random walk assumes that the changes in the parameter are random and without a clear trend, which may not be appropriate in all cases.
#########################################################################
#COMPARACION DE MODELOS con distinto desembarque
#########################################################################
#PLOT labels, name of each model run
legend.labels <- c('herring-aspm',
'herring-just-2019-as-equilibrium')
#read in all model runs
#note if cover=T you need a hessian; if covar=F you do not need a hessian
biglist <- SSgetoutput(keyvec = NULL,
dirvec = c(dir1,dir2,dir3),
getcovar = F)
#create summary of model runs from list above
summaryoutput <- SSsummarize(biglist)
SSplotComparisons(summaryoutput, subplots = 1:20, plot = TRUE,
print = T, endyrvec = "default", indexfleets = NULL,
indexUncertainty = FALSE,
indexQlabel = TRUE,
indexQdigits = 4,
indexSEvec = "default",
indexPlotEach = FALSE,
labels = c("Year", "Spawning biomass (t)",
"Relative spawning biomass",
"Age-0 recruits (1,000s)",
"Recruitment deviations",
"Index",
"Log index",
"1 - SPR",
"Density",
"Management target",
"Minimum stock size threshold",
"Spawning output",
"Harvest rate"), col = NULL,
shadecol = NULL, pch = NULL,
lty = 1, lwd = 2, spacepoints = 10, staggerpoints = 1,
initpoint = 0, tickEndYr = TRUE, shadeForecast = TRUE,
xlim = "default", ylimAdj = 1, xaxs = "r", yaxs = "r",
type = "o", uncertainty = TRUE,
shadealpha = 0.1, legend = TRUE,
legendlabels = "default", legendloc = "topright",
legendorder = "default", legendncol = 1, sprtarg = NULL,
btarg = NULL, minbthresh = NULL,
pwidth = 6.5, pheight = 5,
punits = "in", res = 300, ptsize = 10,
plotdir = "~/DOCAS/Workshops/SS3_ICES_Course",
filenameprefix = "",
densitynames = c("SSB_Virgin", "R0"),
densityxlabs = "default", densityscalex = 1,
densityscaley = 1, densityadjust = 1,
densitysymbols = TRUE,
densitytails = TRUE,
densitymiddle = FALSE,
densitylwd = 1,
fix0 = TRUE, new = TRUE,
add = FALSE,
par = list(mar = c(5, 4, 1, 1) + 0.1),
verbose = TRUE, mcmcVec = FALSE,
show_equilibrium = TRUE)
SStableComparisons(summaryoutput,
likenames = c("TOTAL",
"Survey",
"Length_comp",
"Age_comp",
"priors",
"Size_at_age"),
names = c("Recr_Virgin",
"R0", "steep",
"NatM",
"L_at_Amax",
"VonBert_K",
"SSB_Virg",
"Bratio_2017", "SPRratio_2016"),
digits = NULL, modelnames = "default", csv = FALSE,
csvdir ="~/DOCAS/Workshops/SS3_ICES_Course",
csvfile = "parameter_comparison_table.csv", verbose = TRUE,
mcmc = FALSE)Analisis retrospectivo
#do retrospective model runs
SS_doRetro('~/DOCAS/Workshops/SS3_ICES_Course',
'', newsubdir = "retrospectives",
subdirstart = "retro", years = 0:-5,
overwrite = TRUE, exefile = "ss",
extras = "-nox -nohess",
intern = FALSE, CallType = "system",
RemoveBlocks = FALSE)
retroModels <- SSgetoutput(dirvec=file.path('~/DOCAS/Workshops/SS3_ICES_Course',
"retrospectives",
paste("retro",0:-5,sep="")))
retroSummary <- SSsummarize(retroModels)
endyrvec <- retroSummary$endyrs + 0:-5
SSplotComparisons(retroSummary,
endyrvec=endyrvec,
legendlabels=paste("Data",0:-5,"years"),
plotdir='~/DOCAS/Workshops/SS3_ICES_Course',
plot=TRUE,print=T)
TableCompare <- SStableComparisons(retroSummary,
likenames=like,
names=names,
modelnames=c('B','-1','-2','-3','-4','-5'),
csv=TRUE,
csvfile="RetroRuns.csv",verbose=FALSE)Course Notes
- Lorenzen mortality works just in case you have parametric growth. Used a lot in US in Europe (Cardinale, 2023). About F:
- Both use the same F calculations
- Hybrid best when F is low
- At high F, hybrid’s algorithm to maintain match to catch inhibits movement towards convergence
- Parameter approach can be configured to get starting value for F parameter during early phases using hybrid, then converts to parameter to finish the run
- F method = 4 gives best flexibility because it allows hybrid vs parameter conversion on a fleet-specific basis (today is the best method. If you have fleet specific).
- 0.23 # F ballpark value in units of annual_F. Is usefuul when we have outlier data (High biomas or high F) and you can set what year happen this, i. e. ; -2008 # F ballpark year (neg value to disable).
- Body weigth is in kilograms -Options 3 & 4 are experimental in SS3, but match some other models, like small pelagics stock.
- If fixing without having an intended rationale, then let the estimated recruitments deviate from the curve: –- use recdev option 2, not 1, for more flexibility –- set lambda on recdevs to low value (<0.10)